{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "<center>\n",
    "<img src=\"https://raw.githubusercontent.com/Yorko/mlcourse.ai/master/img/ods_stickers.jpg\">\n",
    "## Open Machine Learning Course\n",
    "<center>\n",
    "Auteur: [Sergey Korolev](https://www.linkedin.com/in/sokorolev/), Ingénieur Logiciel chez Snap Inc. <br>\n",
    "Traduit et édité par [Egor Polusmak](https://www.linkedin.com/in/egor-polusmak/), [Anastasia Manokhina](https://www.linkedin.com/in/anastasiamanokhina/), [Anna Golovchenko](https://www.linkedin.com/in/anna-golovchenko-b0ba5a112/), [Eugene Mashkin](https://www.linkedin.com/in/eugene-mashkin-88490883/), [Yuanyuan Pao](https://www.linkedin.com/in/yuanyuanpao/) et [Ousmane Cissé](https://www.linkedin.com/in/ousmane-cissé/).\n",
    "\n",
    "Ce matériel est soumis aux termes et conditions de la licence [Creative Commons CC BY-NC-SA 4.0](https://creativecommons.org/licenses/by-nc-sa/4.0/). La libre utilisation est autorisée à toute fin non commerciale avec indication obligatoire des noms des auteurs et de la source."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "# <center>Topic 7. Apprentissage non supervisé: PCA et clustering\n",
    "\n",
    "Bienvenue dans la septième partie de notre cours Open Machine Learning!\n",
    "  \n",
    "Dans cette leçon, nous travaillerons avec des méthodes d'apprentissage non supervisées telles que l'analyse en composantes principales (ACP) et le clustering. Vous apprendrez pourquoi et comment nous pouvons réduire la dimensionnalité des données d'origine et quelles sont les principales approches pour regrouper des points de données similaires."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "### Plan de l'article\n",
    "1. Introduction\n",
    "2. PCA\n",
    " - Intuition, théories et problèmes d'application\n",
    " - Exemples d'application\n",
    "3. Analyse des clusters\n",
    " - K-Means\n",
    " - Propagation d'affinité\n",
    " - Regroupement spectral\n",
    " - Cluster agglomératif\n",
    " - Mesures de précision\n",
    "4. Liens utiles"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "## 1. Introduction\n",
    "\n",
    "La principale caractéristique des algorithmes d'apprentissage non supervisés, par rapport aux méthodes de classification et de régression, est que les données d'entrée ne sont pas étiquetées (c'est-à-dire qu'aucune étiquette ou classe n'est donnée) et que l'algorithme apprend la structure des données sans aucune assistance. Cela crée deux différences principales. Tout d'abord, cela nous permet de traiter de grandes quantités de données car les données n'ont pas besoin d'être étiquetées manuellement. Deuxièmement, il est difficile d'évaluer la qualité d'un algorithme non supervisé en raison de l'absence d'une métrique de \"justesse\" explicite telle qu'elle est utilisée dans l'apprentissage supervisé.\n",
    "\n",
    "L'une des tâches les plus courantes dans l'apprentissage non supervisé est la réduction de la dimensionnalité. D'une part, la réduction de la dimensionnalité peut aider à la visualisation des données (par exemple la méthode **t-SNE**) tandis que, d'autre part, elle peut aider à gérer la multicolinéarité de vos données et à préparer les données pour une méthode d'apprentissage supervisé (par exemple, les **arbres de décision**)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "## 2. Analyse en composantes principales (ACP)\n",
    "\n",
    "### Intuition, théories et problèmes d'application\n",
    "\n",
    "L'analyse en composantes principales est l'une des méthodes les plus simples, les plus intuitives et les plus utilisées pour réduire la dimensionnalité, projetant des données sur son sous-espace d'entités orthogonales.\n",
    "\n",
    "\n",
    "<img align=\"right\" src=\"https://habrastorage.org/getpro/habr/post_images/bb6/fe7/f06/bb6fe7f06e114bcc9c354a1cb025b966.png\" width=\"400\">\n",
    "\n",
    "\n",
    "Plus généralement, toutes les observations peuvent être considérées comme un ellipsoïde dans un sous-espace d'un espace caractéristique initial, et la nouvelle base définie dans ce sous-espace est alignée avec les axes ellipsoïdes. Cette hypothèse nous permet de supprimer les caractéristiques hautement corrélées puisque les vecteurs d'ensemble de base sont orthogonaux.\n",
    "Dans le cas général, la dimensionnalité ellipsoïde résultante correspond à la dimensionnalité spatiale initiale, mais l'hypothèse que nos données se trouvent dans un sous-espace avec une dimension plus petite nous permet de couper l'espace \"excessif\" avec la nouvelle projection (sous-espace). Nous accomplissons cela d'une manière «gourmande», en sélectionnant séquentiellement chacun des axes ellipsoïdes en identifiant où la dispersion est maximale.\n",
    " \n",
    "\n",
    "> \"Pour gérer les hyper-plans dans un espace en 14 dimensions, visualisez un espace 3D et dites\" quatorze \"très fort. Tout le monde le fait.\" - Geoffrey Hinton\n",
    "\n",
    "\n",
    "Jetons un coup d'œil à la formulation mathématique de ce processus:\n",
    "\n",
    "Afin de diminuer la dimensionnalité de nos données de $n$ à $k$ avec $k \\leq n$, nous trions notre liste d'axes par ordre de dispersion décroissante et prenons le top-$k$ d'entre eux.\n",
    "\n",
    "Nous commençons par calculer la dispersion et la covariance des caractéristiques initiales. Cela se fait généralement avec la matrice de covariance. Selon la définition de covariance, la covariance de deux entités est calculée comme suit: $$cov(X_i, X_j) = E[(X_i - \\mu_i) (X_j - \\mu_j)] = E[X_i X_j] - \\mu_i \\mu_j,$$ où $\\mu_i$ est la valeur attendue de l'entité $i$th. Il convient de noter que la covariance est symétrique et que la covariance d'un vecteur avec lui-même est égale à sa dispersion.\n",
    "\n",
    "Par conséquent, la matrice de covariance est symétrique avec la dispersion des caractéristiques correspondantes sur la diagonale. Les valeurs non diagonales sont les covariances de la paire d'entités correspondante. En termes de matrices où $\\mathbf{X}$ est la matrice des observations, la matrice de covariance est la suivante:\n",
    "\n",
    "$$\\Sigma = E[(\\mathbf{X} - E[\\mathbf{X}]) (\\mathbf{X} - E[\\mathbf{X}])^{T}]$$\n",
    "\n",
    "Brève récapitulation : les matrices, en tant qu'opérateurs linéaires, ont des valeurs propres et des vecteurs propres. Ils sont très pratiques car ils décrivent des parties de notre espace qui ne tournent pas et ne s'étirent que lorsque nous leur appliquons des opérateurs linéaires; les vecteurs propres restent dans la même direction mais sont étirés par une valeur propre correspondante. Formellement, une matrice $M$ avec vecteur propre $w_i$ et valeur propre $\\lambda_i$ satisfait cette équation: $M w_i = \\lambda_i w_i$.\n",
    "\n",
    "La matrice de covariance d'un échantillon $\\mathbf{X}$ peut être écrite comme un produit de $\\mathbf{X}^{T} \\mathbf{X}$. Selon le [quotient de Rayleigh](https://en.wikipedia.org/wiki/Rayleigh_quotient), la variation maximale de notre échantillon se situe le long du vecteur propre de cette matrice et est cohérente avec la valeur propre maximale. Par conséquent, les principaux composants que nous visons à retenir des données ne sont que les vecteurs propres correspondant aux valeurs propres les plus élevées de la matrice-$k$.\n",
    "\n",
    "Les prochaines étapes sont plus faciles à digérer. Nous multiplions la matrice de nos données $X$ par ces composants pour obtenir la projection de nos données sur la base orthogonale des composants choisis. Si le nombre de composants était inférieur à la dimensionnalité d'espace initiale, n'oubliez pas que nous perdrons certaines informations lors de l'application de cette transformation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "## Exemples\n",
    "### Ensemble de données sur l'iris de Fisher\n",
    "\n",
    "Commençons par télécharger tous les modules essentiels et essayons l'exemple iris de la documentation `scikit-learn`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "import seaborn as sns\n",
    "\n",
    "sns.set(style=\"white\")\n",
    "%matplotlib inline\n",
    "%config InlineBackend.figure_format = 'retina'\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "from sklearn import datasets, decomposition\n",
    "\n",
    "# Loading the dataset\n",
    "iris = datasets.load_iris()\n",
    "X = iris.data\n",
    "y = iris.target\n",
    "\n",
    "# Let's create a beautiful 3d-plot\n",
    "fig = plt.figure(1, figsize=(6, 5))\n",
    "plt.clf()\n",
    "ax = Axes3D(fig, rect=[0, 0, 0.95, 1], elev=48, azim=134)\n",
    "\n",
    "plt.cla()\n",
    "\n",
    "for name, label in [(\"Setosa\", 0), (\"Versicolour\", 1), (\"Virginica\", 2)]:\n",
    "    ax.text3D(\n",
    "        X[y == label, 0].mean(),\n",
    "        X[y == label, 1].mean() + 1.5,\n",
    "        X[y == label, 2].mean(),\n",
    "        name,\n",
    "        horizontalalignment=\"center\",\n",
    "        bbox=dict(alpha=0.5, edgecolor=\"w\", facecolor=\"w\"),\n",
    "    )\n",
    "# Change the order of labels, so that they match\n",
    "y_clr = np.choose(y, [1, 2, 0]).astype(np.float)\n",
    "ax.scatter(X[:, 0], X[:, 1], X[:, 2], c=y_clr, cmap=plt.cm.nipy_spectral)\n",
    "\n",
    "ax.w_xaxis.set_ticklabels([])\n",
    "ax.w_yaxis.set_ticklabels([])\n",
    "ax.w_zaxis.set_ticklabels([]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "Voyons maintenant comment PCA améliorera les résultats d'un modèle simple qui n'est pas capable d'ajuster correctement toutes les données d'entraînement:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.metrics import accuracy_score, roc_auc_score\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.tree import DecisionTreeClassifier\n",
    "\n",
    "# Train, test splits\n",
    "X_train, X_test, y_train, y_test = train_test_split(\n",
    "    X, y, test_size=0.3, stratify=y, random_state=42\n",
    ")\n",
    "\n",
    "# Decision trees with depth = 2\n",
    "clf = DecisionTreeClassifier(max_depth=2, random_state=42)\n",
    "clf.fit(X_train, y_train)\n",
    "preds = clf.predict_proba(X_test)\n",
    "print(\"Accuracy: {:.5f}\".format(accuracy_score(y_test, preds.argmax(axis=1))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "Essayons à nouveau, mais, cette fois, réduisons la dimensionnalité à 2 dimensions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Using PCA from sklearn PCA\n",
    "pca = decomposition.PCA(n_components=2)\n",
    "X_centered = X - X.mean(axis=0)\n",
    "pca.fit(X_centered)\n",
    "X_pca = pca.transform(X_centered)\n",
    "\n",
    "# Plotting the results of PCA\n",
    "plt.plot(X_pca[y == 0, 0], X_pca[y == 0, 1], \"bo\", label=\"Setosa\")\n",
    "plt.plot(X_pca[y == 1, 0], X_pca[y == 1, 1], \"go\", label=\"Versicolour\")\n",
    "plt.plot(X_pca[y == 2, 0], X_pca[y == 2, 1], \"ro\", label=\"Virginica\")\n",
    "plt.legend(loc=0);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Test-train split and apply PCA\n",
    "X_train, X_test, y_train, y_test = train_test_split(\n",
    "    X_pca, y, test_size=0.3, stratify=y, random_state=42\n",
    ")\n",
    "\n",
    "clf = DecisionTreeClassifier(max_depth=2, random_state=42)\n",
    "clf.fit(X_train, y_train)\n",
    "preds = clf.predict_proba(X_test)\n",
    "print(\"Accuracy: {:.5f}\".format(accuracy_score(y_test, preds.argmax(axis=1))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "La précision n'a pas augmenté de manière significative dans ce cas, mais, avec d'autres ensembles de données avec un nombre élevé de dimensions, l'ACP peut considérablement améliorer la précision des arbres de décision et d'autres méthodes d'ensemble."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "Voyons maintenant le pourcentage de variance qui peut être expliqué par chacun des composants sélectionnés."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i, component in enumerate(pca.components_):\n",
    "    print(\n",
    "        \"{} component: {}% of initial variance\".format(\n",
    "            i + 1, round(100 * pca.explained_variance_ratio_[i], 2)\n",
    "        )\n",
    "    )\n",
    "    print(\n",
    "        \" + \".join(\n",
    "            \"%.3f x %s\" % (value, name)\n",
    "            for value, name in zip(component, iris.feature_names)\n",
    "        )\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "### Ensemble de données de nombres manuscrits\n",
    "\n",
    "Regardons l'ensemble de données de nombres manuscrits que nous avons utilisé précédemment dans la [3e leçon](https://habrahabr.ru/company/ods/blog/322534/#derevya-resheniy-i-metod-blizhayshih-sosedey-v-zadache -raspoznavaniya-rukopisnyh-cifr-mnist)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "digits = datasets.load_digits()\n",
    "X = digits.data\n",
    "y = digits.target"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "Commençons par visualiser nos données. Récupérez les 10 premiers chiffres. Les nombres sont représentés par des matrices 8 x 8 avec l'intensité de couleur pour chaque pixel. Chaque matrice est aplatie en un vecteur de 64 nombres, nous obtenons donc la version caractéristique des données."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# f, axes = plt.subplots(5, 2, sharey=True, figsize=(16,6))\n",
    "plt.figure(figsize=(16, 6))\n",
    "for i in range(10):\n",
    "    plt.subplot(2, 5, i + 1)\n",
    "    plt.imshow(X[i, :].reshape([8, 8]), cmap=\"gray\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "Nos données ont 64 dimensions, mais nous allons les réduire à seulement 2 et voir que, même avec seulement 2 dimensions, nous pouvons clairement voir que les chiffres se séparent en grappes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pca = decomposition.PCA(n_components=2)\n",
    "X_reduced = pca.fit_transform(X)\n",
    "\n",
    "print(\"Projecting %d-dimensional data to 2D\" % X.shape[1])\n",
    "\n",
    "plt.figure(figsize=(12, 10))\n",
    "plt.scatter(\n",
    "    X_reduced[:, 0],\n",
    "    X_reduced[:, 1],\n",
    "    c=y,\n",
    "    edgecolor=\"none\",\n",
    "    alpha=0.7,\n",
    "    s=40,\n",
    "    cmap=plt.cm.get_cmap(\"nipy_spectral\", 10),\n",
    ")\n",
    "plt.colorbar()\n",
    "plt.title(\"MNIST. PCA projection\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "En effet, avec t-SNE, l'image semble meilleure car PCA a une contrainte linéaire alors que t-SNE n'en a pas. Cependant, même avec un si petit ensemble de données, l'algorithme t-SNE prend beaucoup plus de temps à compléter que PCA."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "from sklearn.manifold import TSNE\n",
    "\n",
    "tsne = TSNE(random_state=17)\n",
    "\n",
    "X_tsne = tsne.fit_transform(X)\n",
    "\n",
    "plt.figure(figsize=(12, 10))\n",
    "plt.scatter(\n",
    "    X_tsne[:, 0],\n",
    "    X_tsne[:, 1],\n",
    "    c=y,\n",
    "    edgecolor=\"none\",\n",
    "    alpha=0.7,\n",
    "    s=40,\n",
    "    cmap=plt.cm.get_cmap(\"nipy_spectral\", 10),\n",
    ")\n",
    "plt.colorbar()\n",
    "plt.title(\"MNIST. t-SNE projection\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "Dans la pratique, nous choisirions le nombre de composants principaux de façon à pouvoir expliquer 90% de la dispersion initiale des données (via le `explained_variance_ratio`). Ici, cela signifie conserver 21 composants principaux; par conséquent, nous réduisons la dimensionnalité de 64 entités à 21."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pca = decomposition.PCA().fit(X)\n",
    "\n",
    "plt.figure(figsize=(10, 7))\n",
    "plt.plot(np.cumsum(pca.explained_variance_ratio_), color=\"k\", lw=2)\n",
    "plt.xlabel(\"Number of components\")\n",
    "plt.ylabel(\"Total explained variance\")\n",
    "plt.xlim(0, 63)\n",
    "plt.yticks(np.arange(0, 1.1, 0.1))\n",
    "plt.axvline(21, c=\"b\")\n",
    "plt.axhline(0.9, c=\"r\")\n",
    "plt.show();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "## 2. Clustering\n",
    "\n",
    "L'idée principale derrière le clustering est assez simple. Fondamentalement, nous nous disons: «J'ai ces points ici, et je peux voir qu'ils s'organisent en groupes. Ce serait bien de décrire ces choses plus concrètement et, quand un nouveau point arrive, de les affecter au bon groupe . \" Cette idée générale encourage l'exploration et ouvre une variété d'algorithmes pour le clustering.\n",
    "\n",
    "<figure><img align=\"center\" src=\"https://habrastorage.org/getpro/habr/post_images/8b9/ae5/586/8b9ae55861f22a2809e8b3a00ef815ad.png\"><figcaption> *Les exemples des résultats de différents algorithmes de scikit-learn* </figcaption></figure>\n",
    "\n",
    "Les algorithmes répertoriés ci-dessous ne couvrent pas toutes les méthodes de clustering disponibles, mais ce sont les plus couramment utilisées.\n",
    "\n",
    "### K-means\n",
    "\n",
    "L'algorithme K-means est le plus populaire et le plus simple de tous les algorithmes de clustering. Voici comment cela fonctionne:\n",
    "1. Sélectionnez le nombre de grappes $k$ que vous pensez être le nombre optimal.\n",
    "2. Initialisez les points $k$ comme des \"centroïdes\" au hasard dans l'espace de nos données.\n",
    "3. Attribuez chaque observation à son centroïde le plus proche.\n",
    "4. Mettez à jour les centroïdes au centre de l'ensemble des observations attribuées.\n",
    "5. Répétez les étapes 3 et 4 un nombre fixe de fois ou jusqu'à ce que tous les centroïdes soient stables (c'est-à-dire qu'ils ne changent plus à l'étape 4).\n",
    "\n",
    "Cet algorithme est facile à décrire et à visualiser."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let's begin by allocation 3 cluster's points\n",
    "X = np.zeros((150, 2))\n",
    "\n",
    "np.random.seed(seed=42)\n",
    "X[:50, 0] = np.random.normal(loc=0.0, scale=0.3, size=50)\n",
    "X[:50, 1] = np.random.normal(loc=0.0, scale=0.3, size=50)\n",
    "\n",
    "X[50:100, 0] = np.random.normal(loc=2.0, scale=0.5, size=50)\n",
    "X[50:100, 1] = np.random.normal(loc=-1.0, scale=0.2, size=50)\n",
    "\n",
    "X[100:150, 0] = np.random.normal(loc=-1.0, scale=0.2, size=50)\n",
    "X[100:150, 1] = np.random.normal(loc=2.0, scale=0.5, size=50)\n",
    "\n",
    "plt.figure(figsize=(5, 5))\n",
    "plt.plot(X[:, 0], X[:, 1], \"bo\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Scipy has function that takes 2 tuples and return\n",
    "# calculated distance between them\n",
    "from scipy.spatial.distance import cdist\n",
    "\n",
    "# Randomly allocate the 3 centroids\n",
    "np.random.seed(seed=42)\n",
    "centroids = np.random.normal(loc=0.0, scale=1.0, size=6)\n",
    "centroids = centroids.reshape((3, 2))\n",
    "\n",
    "cent_history = []\n",
    "cent_history.append(centroids)\n",
    "\n",
    "for i in range(3):\n",
    "    # Calculating the distance from a point to a centroid\n",
    "    distances = cdist(X, centroids)\n",
    "    # Checking what's the closest centroid for the point\n",
    "    labels = distances.argmin(axis=1)\n",
    "\n",
    "    # Labeling the point according the point's distance\n",
    "    centroids = centroids.copy()\n",
    "    centroids[0, :] = np.mean(X[labels == 0, :], axis=0)\n",
    "    centroids[1, :] = np.mean(X[labels == 1, :], axis=0)\n",
    "    centroids[2, :] = np.mean(X[labels == 2, :], axis=0)\n",
    "\n",
    "    cent_history.append(centroids)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let's plot K-means\n",
    "plt.figure(figsize=(8, 8))\n",
    "for i in range(4):\n",
    "    distances = cdist(X, cent_history[i])\n",
    "    labels = distances.argmin(axis=1)\n",
    "\n",
    "    plt.subplot(2, 2, i + 1)\n",
    "    plt.plot(X[labels == 0, 0], X[labels == 0, 1], \"bo\", label=\"cluster #1\")\n",
    "    plt.plot(X[labels == 1, 0], X[labels == 1, 1], \"co\", label=\"cluster #2\")\n",
    "    plt.plot(X[labels == 2, 0], X[labels == 2, 1], \"mo\", label=\"cluster #3\")\n",
    "    plt.plot(cent_history[i][:, 0], cent_history[i][:, 1], \"rX\")\n",
    "    plt.legend(loc=0)\n",
    "    plt.title(\"Step {:}\".format(i + 1));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "Nous avons utilisé la distance euclidienne, mais l'algorithme convergera avec toute autre métrique. Vous pouvez non seulement faire varier le nombre d'étapes ou les critères de convergence, mais également la mesure de distance entre les points et les centroïdes de cluster.\n",
    "\n",
    "Une autre \"caractéristique\" de cet algorithme est sa sensibilité aux positions initiales des centroïdes du cluster. Vous pouvez exécuter l'algorithme plusieurs fois, puis faire la moyenne de tous les résultats du centroïde."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "## Choix du nombre de clusters pour K-means\n",
    "\n",
    "Contrairement aux tâches d'apprentissage supervisé telles que la classification et la régression, le clustering nécessite plus d'efforts pour choisir le critère d'optimisation. Habituellement, lorsque nous travaillons avec k-means, nous optimisons la somme des distances au carré entre les observations et leurs centroïdes.\n",
    "\n",
    "$$ J(C) = \\sum_{k=1}^K\\sum_{i~\\in~C_k} ||x_i - \\mu_k|| \\rightarrow \\min\\limits_C,$$\n",
    "\n",
    "où $C$ - est un ensemble de clusters avec la puissance $K$, $\\mu_k$ est un centre de gravité d'un cluster $C_k$.\n",
    "\n",
    "Cette définition semble raisonnable - nous voulons que nos observations soient aussi proches que possible de leurs centroïdes. Mais, il y a un problème - l'optimum est atteint lorsque le nombre de centroïdes est égal au nombre d'observations, donc vous vous retrouveriez avec chaque observation comme son propre cluster séparé.\n",
    "\n",
    "Afin d'éviter ce cas, nous devons choisir un certain nombre de clusters après quoi une fonction $J(C_k)$ diminue moins rapidement. Plus formellement,\n",
    "$$ D(k) = \\frac{|J(C_k) - J(C_{k+1})|}{|J(C_{k-1}) - J(C_k)|}  \\rightarrow \\min\\limits_k $$\n",
    "\n",
    "Regardons un exemple."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.cluster import KMeans"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "inertia = []\n",
    "for k in range(1, 8):\n",
    "    kmeans = KMeans(n_clusters=k, random_state=1).fit(X)\n",
    "    inertia.append(np.sqrt(kmeans.inertia_))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(range(1, 8), inertia, marker=\"s\")\n",
    "plt.xlabel(\"$k$\")\n",
    "plt.ylabel(\"$J(C_k)$\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "Nous voyons que $J(C_k)$ diminue de manière significative jusqu'à ce que le nombre de clusters soit de 3 puis ne change plus autant. Cela signifie que le nombre optimal de clusters est de 3."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "#### Problèmes\n",
    "\n",
    "En soi, K-means est NP-difficile. Pour les dimensions $d$, les clusters $k$ et les observations $n$, nous trouverons une solution en temps $O(n^{d k+1})$. Il y a quelques heuristiques pour gérer cela; un exemple est MiniBatch K-means, qui prend des parties (lots) de données au lieu d'ajuster l'ensemble de données, puis déplace les centroïdes en prenant la moyenne des étapes précédentes. Comparez l'implémentation de K-means et MiniBatch K-means dans la [documentation sckit-learn](http://scikit-learn.org/stable/auto_examples/cluster/plot_mini_batch_kmeans.html).\n",
    "\n",
    "L'[implémentation](http://scikit-learn.org/stable/modules/generated/sklearn.cluster.KMeans.html) de l'algorithme utilisant `scikit-learn` a ses avantages tels que la possibilité d'indiquer le nombre d'initialisations avec le paramètre de fonction `n_init`, qui nous permet d'identifier des centroïdes plus robustes. De plus, ces analyses peuvent être effectuées en parallèle pour diminuer le temps de calcul."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "## Propagation d'affinité\n",
    "\n",
    "La propagation d'affinité est un autre exemple d'algorithme de clustering. Contrairement à K-means, cette approche ne nécessite pas de fixer au préalable le nombre de clusters. L'idée principale ici est que nous aimerions regrouper nos données en fonction de la similitude des observations (ou de la façon dont elles «correspondent» les unes aux autres).\n",
    "\n",
    "Définissons une métrique de similitude telle que $s(x_i, x_j) > s(x_i, x_k)$ si une observation $x_i$ est plus similaire à l'observation $x_j$ et moins similaire à l'observation $x_k$. Un exemple simple d'une telle métrique de similitude est un carré négatif de la distance $s(x_i, x_j) = - ||x_i - x_j||^{2}$.\n",
    "\n",
    "\n",
    "Décrivons maintenant la \"correspondance\" en faisant deux matrices nulles. L'un d'eux, $r_{i,k}$, détermine dans quelle mesure l'observation $k$th est en tant que «modèle de rôle» pour l'observation $i$th par rapport à tous les autres «modèles de rôle» possibles. Une autre matrice, $a_{i,k}$, détermine dans quelle mesure il serait approprié pour $i$th observation de prendre l'observation $k$th comme \"modèle\".\n",
    "\n",
    "Les matrices sont mises à jour séquentiellement avec les règles suivantes:\n",
    "\n",
    "$$r_{i,k} \\leftarrow s_(x_i, x_k) - \\max_{k' \\neq k} \\left\\{ a_{i,k'} + s(x_i, x_k') \\right\\}$$\n",
    "\n",
    "$$a_{i,k} \\leftarrow \\min \\left( 0, r_{k,k} + \\sum_{i' \\not\\in \\{i,k\\}} \\max(0, r_{i',k}) \\right), \\ \\ \\  i \\neq k$$\n",
    "\n",
    "$$a_{k,k} \\leftarrow \\sum_{i' \\neq k} \\max(0, r_{i',k})$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "## Regroupement spectral\n",
    "\n",
    "Le regroupement spectral combine certaines des approches décrites ci-dessus pour créer une méthode de regroupement plus solide.\n",
    "\n",
    "Tout d'abord, cet algorithme nous oblige à définir la matrice de similitude pour les observations appelée matrice d'adjacence. Cela peut être fait d'une manière similaire à celle de l'algorithme de propagation d'affinité: $A_{i, j} = - ||x_i - x_j||^{2}$. Cette matrice décrit un graphique complet avec les observations sous forme de sommets et la valeur de similitude estimée entre une paire d'observations sous forme de poids de bord pour cette paire de sommets. Pour la métrique définie ci-dessus et les observations bidimensionnelles, cela est assez intuitif - deux observations sont similaires si le bord entre elles est plus court.\n",
    "Nous aimerions diviser le graphique en deux sous-graphiques de telle sorte que chaque observation de chaque sous-graphique soit similaire à une autre observation de ce sous-graphique. Formellement, il s'agit d'un problème de coupes normalisées; pour plus de détails, nous vous recommandons de lire [ce document](http://people.eecs.berkeley.edu/~malik/papers/SM-ncut.pdf)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "## Cluster agglomératif"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "L'algorithme suivant est le plus simple et le plus facile à comprendre parmi tous les algorithmes de clustering sans nombre fixe de clusters.\n",
    "\n",
    "\n",
    "L'algorithme est assez simple:\n",
    "1. Nous commençons par assigner chaque observation à son propre cluster\n",
    "2. Ensuite, triez les distances par paires entre les centres des grappes dans l'ordre décroissant\n",
    "3. Prenez les deux grappes voisines les plus proches et fusionnez-les, puis recalculez les centres\n",
    "4. Répétez les étapes 2 et 3 jusqu'à ce que toutes les données soient fusionnées en un seul cluster\n",
    "\n",
    "Le processus de recherche du cluster le plus proche peut être mené avec différentes méthodes de délimitation des observations:\n",
    "1. Liaison unique\n",
    "$d(C_i, C_j) = min_{x_i \\in C_i, x_j \\in C_j} ||x_i - x_j||$\n",
    "2. Lien complet\n",
    "$d(C_i, C_j) = max_{x_i \\in C_i, x_j \\in C_j} ||x_i - x_j||$\n",
    "3. Lien moyen\n",
    "$d(C_i, C_j) = \\frac{1}{n_i n_j} \\sum_{x_i \\in C_i} \\sum_{x_j \\in C_j} ||x_i - x_j||$\n",
    "4. Liaison centroïde\n",
    "$d(C_i, C_j) = ||\\mu_i - \\mu_j||$\n",
    "\n",
    "Le 3ème est le plus efficace en temps de calcul car il ne nécessite pas de recalculer les distances à chaque fusion des clusters.\n",
    "\n",
    "Les résultats peuvent être visualisés comme un bel arbre de cluster (dendogramme) pour aider à reconnaître le moment où l'algorithme doit être arrêté pour obtenir des résultats optimaux. Il existe de nombreux outils Python pour construire ces dendogrammes pour le clustering agglomératif.\n",
    "\n",
    "Prenons un exemple avec les clusters que nous avons obtenus de K-means:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.cluster import hierarchy\n",
    "from scipy.spatial.distance import pdist\n",
    "\n",
    "X = np.zeros((150, 2))\n",
    "\n",
    "np.random.seed(seed=42)\n",
    "X[:50, 0] = np.random.normal(loc=0.0, scale=0.3, size=50)\n",
    "X[:50, 1] = np.random.normal(loc=0.0, scale=0.3, size=50)\n",
    "\n",
    "X[50:100, 0] = np.random.normal(loc=2.0, scale=0.5, size=50)\n",
    "X[50:100, 1] = np.random.normal(loc=-1.0, scale=0.2, size=50)\n",
    "\n",
    "X[100:150, 0] = np.random.normal(loc=-1.0, scale=0.2, size=50)\n",
    "X[100:150, 1] = np.random.normal(loc=2.0, scale=0.5, size=50)\n",
    "\n",
    "# pdist will calculate the upper triangle of the pairwise distance matrix\n",
    "distance_mat = pdist(X)\n",
    "# linkage — is an implementation of agglomerative algorithm\n",
    "Z = hierarchy.linkage(distance_mat, \"single\")\n",
    "plt.figure(figsize=(10, 5))\n",
    "dn = hierarchy.dendrogram(Z, color_threshold=0.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "## Mesures de précision\n",
    "\n",
    "Contrairement à la classification, il est difficile d'évaluer la qualité des résultats du clustering. Ici, une métrique ne peut pas dépendre des étiquettes mais uniquement de la qualité de la division. Deuxièmement, nous n'avons généralement pas de véritables étiquettes des observations lorsque nous utilisons le clustering.\n",
    "\n",
    "Il existe des mesures de qualité *internes* et *externes*. Les métriques externes utilisent les informations sur la véritable séparation connue tandis que les métriques internes n'utilisent aucune information externe et évaluent la qualité des clusters en se basant uniquement sur les données initiales. Le nombre optimal de clusters est généralement défini par rapport à certaines métriques internes.\n",
    "\n",
    "Toutes les métriques décrites ci-dessous sont implémentées dans `sklearn.metrics`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "**Adjusted Rand Index (ARI)**\n",
    "\n",
    "Ici, nous supposons que les véritables étiquettes des objets sont connues. Cette métrique ne dépend pas des valeurs des étiquettes mais de la division du cluster de données. Soit $N$ le nombre d'observations dans un échantillon. Soit $a$ le nombre de paires d'observation avec les mêmes étiquettes et situées dans le même cluster, et que $b$ soit le nombre de paires d'observation avec des étiquettes différentes et situées dans différents clusters. L'indice Rand peut être calculé à l'aide de la formule suivante: $$\\text{RI} = \\frac{2(a + b)}{n(n-1)}.$$\n",
    "En d'autres termes, il évalue une part de paires d'observation pour lesquelles ces divisions (résultat initial et clustering) sont cohérentes. L'indice Rand (RI) évalue la similitude des deux divisions du même échantillon. Pour que cet indice soit proche de zéro pour tout résultat de clustering avec n'importe quel $n$ et nombre de clusters, il est essentiel de le mettre à l'échelle, d'où l'indice Adjusted Rand Index: $$\\text{ARI} = \\frac{\\text{RI} - E[\\text{RI}]}{\\max(\\text{RI}) - E[\\text{RI}]}.$$\n",
    "\n",
    "Cette métrique est symétrique et ne dépend pas de la permutation d'étiquette. Par conséquent, cet indice est une mesure des distances entre différentes divisions d'échantillon. $\\text{ARI}$ prend des valeurs dans la plage $[-1, 1]$. Les valeurs négatives indiquent l'indépendance des divisions et les valeurs positives indiquent que ces divisions sont cohérentes (elles correspondent à $\\text{ARI} = 1$)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "**Adjusted Mutual Information (AMI)**\n",
    "\n",
    "Cette métrique est similaire à $\\text{ARI}$. Elle est également symétrique et ne dépend pas des valeurs et de la permutation des étiquettes. Il est défini par la fonction [entropie](https://en.wikipedia.org/wiki/Entropy_(information_theory) et interprète un échantillon divisé comme une distribution discrète (la probabilité d'affecter à un cluster est égal au pourcentage d'objets qu'il contient.) L'index $MI$ est défini comme l'[information mutuelle](https://en.wikipedia.org/wiki/Mutual_information) pour deux distributions, correspondant à l'échantillon divisé en grappes. Intuitivement, les informations mutuelles mesurent la part des informations communes aux deux clusters, c'est-à-dire comment les informations sur l'un d'entre eux diminuent l'incertitude de l'autre.\n",
    "\n",
    "De la même manière que $\\text{ARI}$, $\\text{AMI}$ est défini. Cela nous permet de nous débarrasser de l'augmentation de l'indice $MI$ avec le nombre de clusters. Le $\\text{AMI}$ se situe dans la gamme $[0, 1]$. Des valeurs proches de zéro signifient que les divisions sont indépendantes, et celles proches de 1 signifient qu'elles sont similaires (avec correspondance complète à $\\text{AMI} = 1$)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "**Homogénéité, exhaustivité, V-mesure**\n",
    "\n",
    "Formellement, ces métriques sont également définies en fonction de la fonction d'entropie et de la fonction d'entropie conditionnelle, interprétant les répartitions d'échantillon comme des distributions discrètes: $$h = 1 - \\frac{H(C\\mid K)}{H(C)}, c = 1 - \\frac{H(K\\mid C)}{H(K)},$$\n",
    "où $K$ est un résultat de clustering et $C$ est la division initiale. Par conséquent, $h$ évalue si chaque cluster est composé des mêmes objets de classe, et $c$ mesure dans quelle mesure les mêmes objets de classe s'adaptent aux clusters. Ces métriques ne sont pas symétriques. Les deux se situent dans la plage $[0, 1]$, et des valeurs plus proches de 1 indiquent des résultats de clustering plus précis. Les valeurs de ces métriques ne sont pas mises à l'échelle comme le sont les métriques $\\text{ARI}$ ou $\\text{AMI}$ et dépendent donc du nombre de clusters. Un résultat de clustering aléatoire n'aura pas de valeurs de mesures plus proches de zéro lorsque le nombre de clusters est suffisamment grand et le nombre d'objets est petit. Dans un tel cas, il serait plus raisonnable d'utiliser $\\text{ARI}$. Cependant, avec un grand nombre d'observations (plus de 100) et un nombre de grappes inférieur à 10, ce problème est moins critique et peut être ignoré.\n",
    "\n",
    "$V$-mesure (𝑉-measure) est une combinaison de $h$ et $c$ et est leur moyenne harmonique:\n",
    "$$v = 2\\frac{hc}{h+c}.$$\n",
    "Il est symétrique et mesure la cohérence de deux résultats de regroupement."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "**Silhouette**\n",
    "\n",
    "Contrairement aux métriques décrites ci-dessus, ce coefficient n'implique pas la connaissance des véritables étiquettes des objets. Il nous permet d'estimer la qualité du clustering en utilisant uniquement l'échantillon initial sans étiquette et le résultat du clustering. Pour commencer, pour chaque observation, le coefficient de silhouette est calculé. Soit $a$ la moyenne de la distance entre un objet et d'autres objets au sein d'un cluster et $b$ la distance moyenne d'un objet aux objets du cluster le plus proche (différente de celle à laquelle appartient l'objet). Ensuite, la mesure de silhouette pour cet objet est $$s = \\frac{b - a}{\\max(a, b)}.$$\n",
    "\n",
    "La silhouette d'un échantillon est une valeur moyenne des valeurs de silhouette de cet échantillon. Par conséquent, la distance de la silhouette montre dans quelle mesure la distance entre les objets d'une même classe diffère de la distance moyenne entre les objets de différents groupes. Ce coefficient prend des valeurs dans la plage $[-1, 1]$. Des valeurs proches de -1 correspondent à de mauvais résultats de clustering tandis que des valeurs plus proches de 1 correspondent à des clusters denses et bien définis. Par conséquent, plus la valeur de silhouette est élevée, meilleurs sont les résultats du clustering.\n",
    "\n",
    "À l'aide de silhouette, nous pouvons identifier le nombre optimal de grappes $k$ (si nous ne le savons pas déjà à partir des données) en prenant le nombre de grappes qui maximise le coefficient de silhouette."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "Pour conclure, examinons comment ces mesures fonctionnent avec le jeu de données de nombres manuscrits MNIST:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "from sklearn import datasets, metrics\n",
    "from sklearn.cluster import (\n",
    "    AffinityPropagation,\n",
    "    AgglomerativeClustering,\n",
    "    KMeans,\n",
    "    SpectralClustering,\n",
    ")\n",
    "\n",
    "data = datasets.load_digits()\n",
    "X, y = data.data, data.target\n",
    "\n",
    "algorithms = []\n",
    "algorithms.append(KMeans(n_clusters=10, random_state=1))\n",
    "algorithms.append(AffinityPropagation())\n",
    "algorithms.append(\n",
    "    SpectralClustering(n_clusters=10, random_state=1, affinity=\"nearest_neighbors\")\n",
    ")\n",
    "algorithms.append(AgglomerativeClustering(n_clusters=10))\n",
    "\n",
    "data = []\n",
    "for algo in algorithms:\n",
    "    algo.fit(X)\n",
    "    data.append(\n",
    "        (\n",
    "            {\n",
    "                \"ARI\": metrics.adjusted_rand_score(y, algo.labels_),\n",
    "                \"AMI\": metrics.adjusted_mutual_info_score(y, algo.labels_),\n",
    "                \"Homogenity\": metrics.homogeneity_score(y, algo.labels_),\n",
    "                \"Completeness\": metrics.completeness_score(y, algo.labels_),\n",
    "                \"V-measure\": metrics.v_measure_score(y, algo.labels_),\n",
    "                \"Silhouette\": metrics.silhouette_score(X, algo.labels_),\n",
    "            }\n",
    "        )\n",
    "    )\n",
    "\n",
    "results = pd.DataFrame(\n",
    "    data=data,\n",
    "    columns=[\"ARI\", \"AMI\", \"Homogenity\", \"Completeness\", \"V-measure\", \"Silhouette\"],\n",
    "    index=[\"K-means\", \"Affinity\", \"Spectral\", \"Agglomerative\"],\n",
    ")\n",
    "\n",
    "results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "## 4. Liens utiles\n",
    "- Overview of clustering methods in the [scikit-learn doc](http://scikit-learn.org/stable/modules/clustering.html).\n",
    "- [Q&A](http://stats.stackexchange.com/questions/2691/making-sense-of-principal-component-analysis-eigenvectors-eigenvalues) for PCA with examples\n",
    "- [Notebook](https://github.com/diefimov/MTH594_MachineLearning/blob/master/ipython/Lecture10.ipynb) on k-means and the EM-algorithm</div><i class=\"fa fa-lightbulb-o \"></i>"
   ]
  }
 ],
 "metadata": {
  "hide_input": false,
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.4"
  },
  "nbTranslate": {
   "displayLangs": [
    "*"
   ],
   "hotkey": "alt-t",
   "langInMainMenu": true,
   "sourceLang": "en",
   "targetLang": "fr",
   "useGoogleTranslate": true
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": false,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
